///////////////////////////////////////////////////////////////////////////////
//      Copyright Christopher Kormanyos 2012 - 2015, 2020.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)
//

// This example uses Boost.Multiprecision to implement
// a high-precision Gauss-Laguerre quadrature integration.
// The quadrature is used to calculate the airy_ai(x) function
// for real-valued x on the positive axis with x.ge.1.

// In this way, the integral representation could be seen
// as part of a scheme to calculate real-valued Airy functions
// on the positive axis for medium to large argument.
// A Taylor series or hypergeometric function (not part
// of this example) could be used for smaller arguments.

// This example has been tested with decimal digits counts
// ranging from 21...301, by adjusting the parameter
// local::my_digits10 at compile time.

// The quadrature integral representaion of airy_ai(x) used
// in this example can be found in:

// A. Gil, J. Segura, N.M. Temme, "Numerical Methods for Special
// Functions" (SIAM Press 2007), Sect. 5.3.3, in particular Eq. 5.110,
// page 145. Subsequently, Gil et al's book cites the another work:
// W. Gautschi, "Computation of Bessel and Airy functions and of
// related Gaussian quadrature formulae", BIT, 42 (2002), pp. 110-118.

#include <cmath>
#include <cstdint>
#include <functional>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <sstream>
#include <tuple>
#include <vector>

#include <boost/cstdfloat.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/roots.hpp>
#include <nil/crypto3/multiprecision/cpp_dec_float.hpp>
#include <boost/noncopyable.hpp>

namespace gauss {
    namespace laguerre {

        namespace util {

            void progress_bar(std::ostream& os, const float percent);

            void progress_bar(std::ostream& os, const float percent) {
                std::stringstream strstrm;

                strstrm.precision(1);

                strstrm << std::fixed << percent << "%";

                os << strstrm.str() << "\r";

                os.flush();
            }

        }    // namespace util

        namespace detail {
            template<typename T>
            class laguerre_l_object BOOST_FINAL {
            public:
                laguerre_l_object(const int n, const T a) BOOST_NOEXCEPT : order(n), alpha(a), p1(0), d2(0) {
                }

                laguerre_l_object& operator=(const laguerre_l_object& other) {
                    if (this != other) {
                        order = other.order;
                        alpha = other.alpha;
                        p1 = other.p1;
                        d2 = other.d2;
                    }

                    return *this;
                }

                T operator()(const T& x) const BOOST_NOEXCEPT {
                    // Calculate (via forward recursion):
                    // * the value of the Laguerre function L(n, alpha, x), called (p2),
                    // * the value of the derivative of the Laguerre function (d2),
                    // * and the value of the corresponding Laguerre function of
                    //   previous order (p1).

                    // Return the value of the function (p2) in order to be used as a
                    // function object with Boost.Math root-finding. Store the values
                    // of the Laguerre function derivative (d2) and the Laguerre function
                    // of previous order (p1) in class members for later use.

                    p1 = T(0);
                    T p2 = T(1);
                    d2 = T(0);

                    T j_plus_alpha = alpha;
                    T two_j_plus_one_plus_alpha_minus_x = (1 + alpha) - x;

                    const T my_two = 2;

                    for (int j = 0; j < order; ++j) {
                        const T p0(p1);

                        // Set the value of the previous Laguerre function.
                        p1 = p2;

                        // Use a recurrence relation to compute the value of the Laguerre function.
                        p2 = ((two_j_plus_one_plus_alpha_minus_x * p1) - (j_plus_alpha * p0)) / (j + 1);

                        ++j_plus_alpha;
                        two_j_plus_one_plus_alpha_minus_x += my_two;
                    }

                    // Set the value of the derivative of the Laguerre function.
                    d2 = ((p2 * order) - (j_plus_alpha * p1)) / x;

                    // Return the value of the Laguerre function.
                    return p2;
                }

                const T previous() const BOOST_NOEXCEPT {
                    return p1;
                }
                const T derivative() const BOOST_NOEXCEPT {
                    return d2;
                }

                static bool root_tolerance(const T& a, const T& b) BOOST_NOEXCEPT {
                    using std::fabs;

                    // The relative tolerance here is: ((a - b) * 2) / (a + b).
                    return ((fabs(a - b) * 2) < (fabs(a + b) * std::numeric_limits<T>::epsilon()));
                }

            private:
                const int order;
                const T alpha;
                mutable T p1;
                mutable T d2;
            };

            template<typename T>
            class abscissas_and_weights : private boost::noncopyable {
            public:
                abscissas_and_weights(const int n, const T a) : order(n), alpha(a), xi(), wi() {
                    if (alpha < -20.0F) {
                        // Users can decide to perform different range checking.
                        std::cout << "Range error: the order of the Laguerre function must exceed -20.0." << std::endl;
                    } else {
                        calculate();
                    }
                }

                const std::vector<T>& abscissa_n() const BOOST_NOEXCEPT {
                    return xi;
                }
                const std::vector<T>& weight_n() const BOOST_NOEXCEPT {
                    return wi;
                }

            private:
                const int order;
                const T alpha;

                std::vector<T> xi;
                std::vector<T> wi;

                abscissas_and_weights() : order(), alpha(), xi(), wi() {
                }

                void calculate() {
                    using std::abs;

                    std::cout << "Finding the approximate roots..." << std::endl;

                    std::vector<std::tuple<T, T>> root_estimates;

                    root_estimates.reserve(static_cast<typename std::vector<std::tuple<T, T>>::size_type>(order));

                    const laguerre_l_object<T> laguerre_root_object(order, alpha);

                    // Set the initial values of the step size and the running step
                    // to be used for finding the estimate of the first root.
                    T step_size = 0.01F;
                    T step = step_size;

                    T first_laguerre_root = 0.0F;

                    if (alpha < -1.0F) {
                        // Iteratively step through the Laguerre function using a
                        // small step-size in order to find a rough estimate of
                        // the first zero.

                        const bool this_laguerre_value_is_negative = (laguerre_root_object(T(0)) < 0);

                        BOOST_CONSTEXPR_OR_CONST int j_max = 10000;

                        int j = 0;

                        while ((j < j_max) && (this_laguerre_value_is_negative != (laguerre_root_object(step) < 0))) {
                            // Increment the step size until the sign of the Laguerre function
                            // switches. This indicates a zero-crossing, signalling the next root.
                            step += step_size;

                            ++j;
                        }
                    } else {
                        // Calculate an estimate of the 1st root of a generalized Laguerre
                        // function using either a Taylor series or an expansion in Bessel
                        // function zeros. The Bessel function zeros expansion is from Tricomi.

                        // Here, we obtain an estimate of the first zero of cyl_bessel_j(alpha, x).

                        T j_alpha_m1;

                        if (alpha < 1.4F) {
                            // For small alpha, use a short series obtained from Mathematica(R).
                            // Series[BesselJZero[v, 1], {v, 0, 3}]
                            // N[%, 12]
                            j_alpha_m1 = (((T(0.09748661784476F) * alpha - T(0.17549359276115F)) * alpha +
                                           T(1.54288974259931F)) *
                                              alpha +
                                          T(2.40482555769577F));
                        } else {
                            // For larger alpha, use the first line of Eqs. 10.21.40 in the NIST Handbook.
                            const T alpha_pow_third(boost::math::cbrt(alpha));
                            const T alpha_pow_minus_two_thirds(T(1) / (alpha_pow_third * alpha_pow_third));

                            j_alpha_m1 = alpha * (((((+T(0.043F) * alpha_pow_minus_two_thirds - T(0.0908F)) *
                                                         alpha_pow_minus_two_thirds -
                                                     T(0.00397F)) *
                                                        alpha_pow_minus_two_thirds +
                                                    T(1.033150F)) *
                                                       alpha_pow_minus_two_thirds +
                                                   T(1.8557571F)) *
                                                      alpha_pow_minus_two_thirds +
                                                  T(1.0F));
                        }

                        const T vf = (T(order * 4U) + T(alpha * 2U) + T(2U));
                        const T vf2 = vf * vf;
                        const T j_alpha_m1_sqr = j_alpha_m1 * j_alpha_m1;

                        first_laguerre_root =
                            (j_alpha_m1_sqr * (-T(0.6666666666667F) + ((T(0.6666666666667F) * alpha) * alpha) +
                                               (T(0.3333333333333F) * j_alpha_m1_sqr) + vf2)) /
                            (vf2 * vf);
                    }

                    bool this_laguerre_value_is_negative = (laguerre_root_object(T(0)) < 0);

                    // Re-set the initial value of the step-size based on the
                    // estimate of the first root.
                    step_size = first_laguerre_root / 2;
                    step = step_size;

                    // Step through the Laguerre function using a step-size
                    // of dynamic width in order to find the zero crossings
                    // of the Laguerre function, providing rough estimates
                    // of the roots. Refine the brackets with a few bisection
                    // steps, and store the results as bracketed root estimates.

                    while (root_estimates.size() < static_cast<std::size_t>(order)) {
                        // Increment the step size until the sign of the Laguerre function
                        // switches. This indicates a zero-crossing, signalling the next root.
                        step += step_size;

                        if (this_laguerre_value_is_negative != (laguerre_root_object(step) < 0)) {
                            // We have found the next zero-crossing.

                            // Change the running sign of the Laguerre function.
                            this_laguerre_value_is_negative = (!this_laguerre_value_is_negative);

                            // We have found the first zero-crossing. Put a loose bracket around
                            // the root using a window. Here, we know that the first root lies
                            // between (x - step_size) < root < x.

                            // Before storing the approximate root, perform a couple of
                            // bisection steps in order to tighten up the root bracket.
                            boost::uintmax_t a_couple_of_iterations = 4U;

                            const std::pair<T, T> root_estimate_bracket =
                                boost::math::tools::bisect(laguerre_root_object,
                                                           step - step_size,
                                                           step,
                                                           laguerre_l_object<T>::root_tolerance,
                                                           a_couple_of_iterations);

                            static_cast<void>(a_couple_of_iterations);

                            // Store the refined root estimate as a bracketed range in a tuple.
                            root_estimates.push_back(std::tuple<T, T>(std::get<0>(root_estimate_bracket),
                                                                      std::get<1>(root_estimate_bracket)));

                            if ((root_estimates.size() == 1U) || ((root_estimates.size() % 8U) == 0U) ||
                                (root_estimates.size() == static_cast<std::size_t>(order))) {
                                const float progress =
                                    (100.0F * static_cast<float>(root_estimates.size())) / static_cast<float>(order);

                                std::cout << "root_estimates.size(): " << root_estimates.size() << ", ";

                                util::progress_bar(std::cout, progress);
                            }

                            if (root_estimates.size() >= static_cast<std::size_t>(2U)) {
                                // Determine the next step size. This is based on the distance between
                                // the previous two roots, whereby the estimates of the previous roots
                                // are computed by taking the average of the lower and upper range of
                                // the root-estimate bracket.

                                const T r0 = (std::get<0>(*(root_estimates.crbegin() + 1U)) +
                                              std::get<1>(*(root_estimates.crbegin() + 1U))) /
                                             2;

                                const T r1 =
                                    (std::get<0>(*root_estimates.crbegin()) + std::get<1>(*root_estimates.crbegin())) /
                                    2;

                                const T distance_between_previous_roots = r1 - r0;

                                step_size = distance_between_previous_roots / 3;
                            }
                        }
                    }

                    const T norm_g =
                        ((alpha == 0) ? T(-1) :
                                        -boost::math::tgamma(alpha + order) / boost::math::factorial<T>(order - 1));

                    xi.reserve(root_estimates.size());
                    wi.reserve(root_estimates.size());

                    std::cout << std::endl;

                    // Calculate the abscissas and weights to full precision.
                    for (std::size_t i = static_cast<std::size_t>(0U); i < root_estimates.size(); ++i) {
                        if (((i % 8U) == 0U) || (i == root_estimates.size() - 1U)) {
                            const float progress =
                                (100.0F * static_cast<float>(i + 1U)) / static_cast<float>(root_estimates.size());

                            std::cout << "Calculating abscissas and weights. Processed " << (i + 1U) << ", ";

                            util::progress_bar(std::cout, progress);
                        }

                        // Calculate the abscissas using iterative root-finding.

                        // Select the maximum allowed iterations, being at least 20.
                        // The determination of the maximum allowed iterations is
                        // based on the number of decimal digits in the numerical
                        // type T.
                        BOOST_CONSTEXPR_OR_CONST int local_math_tools_digits10 =
                            static_cast<int>(static_cast<boost::float_least32_t>(boost::math::tools::digits<T>()) *
                                             BOOST_FLOAT32_C(0.301));

                        const boost::uintmax_t number_of_iterations_allowed =
                            (std::max)(20, local_math_tools_digits10 / 2);

                        boost::uintmax_t number_of_iterations_used = number_of_iterations_allowed;

                        // Perform the root-finding using ACM TOMS 748 from Boost.Math.
                        const std::pair<T, T> laguerre_root_bracket =
                            boost::math::tools::toms748_solve(laguerre_root_object,
                                                              std::get<0>(root_estimates.at(i)),
                                                              std::get<1>(root_estimates.at(i)),
                                                              laguerre_l_object<T>::root_tolerance,
                                                              number_of_iterations_used);

                        static_cast<void>(number_of_iterations_used);

                        // Compute the Laguerre root as the average of the values from
                        // the solved root bracket.
                        const T laguerre_root =
                            (std::get<0>(laguerre_root_bracket) + std::get<1>(laguerre_root_bracket)) / 2;

                        // Calculate the weight for this Laguerre root. Here, we calculate
                        // the derivative of the Laguerre function and the value of the
                        // previous Laguerre function on the x-axis at the value of this
                        // Laguerre root.
                        static_cast<void>(laguerre_root_object(laguerre_root));

                        // Store the abscissa and weight for this index.
                        xi.push_back(laguerre_root);
                        wi.push_back(norm_g /
                                     ((laguerre_root_object.derivative() * order) * laguerre_root_object.previous()));
                    }

                    std::cout << std::endl;
                }
            };

            template<typename T>
            struct airy_ai_object BOOST_FINAL {
            public:
                airy_ai_object(const T& x) :
                    my_x(x), my_zeta(((sqrt(x) * x) * 2) / 3), my_factor(make_factor(my_zeta)) {
                }

                T operator()(const T& t) const {
                    using std::sqrt;

                    return my_factor / sqrt(boost::math::cbrt(2 + (t / my_zeta)));
                }

            private:
                const T my_x;
                const T my_zeta;
                const T my_factor;

                airy_ai_object() : my_x(), my_zeta(), my_factor() {
                }

                static T make_factor(const T& z) {
                    using std::exp;
                    using std::sqrt;

                    return 1 / ((sqrt(boost::math::constants::pi<T>()) * sqrt(boost::math::cbrt(z * 48))) *
                                (exp(z) * boost::math::tgamma(T(5) / 6)));
                }
            };
        }    // namespace detail

    }    // namespace laguerre
}    // namespace gauss

// A float_type is created to handle the desired number of decimal digits from `cpp_dec_float` without using
// __expression_templates.
struct local {
    BOOST_STATIC_CONSTEXPR unsigned int my_digits10 = 101U;

    typedef nil::crypto3::multiprecision::number<nil::crypto3::multiprecision::cpp_dec_float<my_digits10>,
                                                 nil::crypto3::multiprecision::et_off>
        float_type;
};

BOOST_STATIC_ASSERT_MSG(local::my_digits10 > 20U,
                        "Error: This example is intended to have more than 20 decimal digits");

int main() {
    // Use Gauss-Laguerre quadrature integration to compute airy_ai(x / 7)
    // with 7 <= x <= 120 and where x is incremented in steps of 1.

    // During development of this example, we have empirically found
    // the numbers of Gauss-Laguerre coefficients needed for convergence
    // when using various counts of base-10 digits.

    // Let's calibrate, for instance, the number of coefficients needed
    // at the point x = 1.

    // Empirical data lead to:
    // Fit[{{21.0, 3.5}, {51.0, 11.1}, {101.0, 22.5}, {201.0, 46.8}}, {1, d, d^2}, d]
    // FullSimplify[%]
    // -1.28301 + (0.235487 + 0.0000178915 d) d

    // We need significantly more coefficients at smaller real values than are needed
    // at larger real values because the slope derivative of airy_ai(x) gets more
    // steep as x approaches zero.

    // This Gauss-Laguerre quadrature is designed for airy_ai(x) with real-valued x >= 1.

    BOOST_CONSTEXPR_OR_CONST boost::float_least32_t d =
        static_cast<boost::float_least32_t>(std::numeric_limits<local::float_type>::digits10);

    BOOST_CONSTEXPR_OR_CONST boost::float_least32_t laguerre_order_factor =
        -1.28301F + ((0.235487F + (0.0000178915F * d)) * d);

    BOOST_CONSTEXPR_OR_CONST int laguerre_order = static_cast<int>(laguerre_order_factor * d);

    std::cout << "std::numeric_limits<local::float_type>::digits10: "
              << std::numeric_limits<local::float_type>::digits10 << std::endl;
    std::cout << "laguerre_order: " << laguerre_order << std::endl;

    typedef gauss::laguerre::detail::abscissas_and_weights<local::float_type> abscissas_and_weights_type;

    const abscissas_and_weights_type the_abscissas_and_weights(laguerre_order, local::float_type(-1) / 6);

    bool result_is_ok = true;

    for (std::uint32_t u = 7U; u < 121U; ++u) {
        const local::float_type x = local::float_type(u) / 7;

        typedef gauss::laguerre::detail::airy_ai_object<local::float_type> airy_ai_object_type;

        const airy_ai_object_type the_airy_ai_object(x);

        const local::float_type airy_ai_value =
            std::inner_product(the_abscissas_and_weights.abscissa_n().cbegin(),
                               the_abscissas_and_weights.abscissa_n().cend(),
                               the_abscissas_and_weights.weight_n().cbegin(),
                               local::float_type(0U),
                               std::plus<local::float_type>(),
                               [&the_airy_ai_object](const local::float_type& this_abscissa,
                                                     const local::float_type& this_weight) -> local::float_type {
                                   return the_airy_ai_object(this_abscissa) * this_weight;
                               });

        static const local::float_type one_third = 1.0F / local::float_type(3U);

        static const local::float_type one_over_pi_times_one_over_sqrt_three =
            sqrt(one_third) * boost::math::constants::one_div_pi<local::float_type>();

        const local::float_type sqrt_x = sqrt(x);

        const local::float_type airy_ai_control =
            (sqrt_x * one_over_pi_times_one_over_sqrt_three) *
            boost::math::cyl_bessel_k(one_third, ((2.0F * x) * sqrt_x) * one_third);

        std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10)
                  << "airy_ai_value  : " << airy_ai_value << std::endl;

        std::cout << std::setprecision(std::numeric_limits<local::float_type>::digits10)
                  << "airy_ai_control: " << airy_ai_control << std::endl;

        const local::float_type delta = fabs(1.0F - (airy_ai_control / airy_ai_value));

        static const local::float_type tol(
            "1E-" + boost::lexical_cast<std::string>(std::numeric_limits<local::float_type>::digits10 - 7U));

        result_is_ok &= (delta < tol);
    }

    std::cout << std::endl << "Total... result_is_ok: " << std::boolalpha << result_is_ok << std::endl;
}    // int main()

/*


Partial output:

//[gauss_laguerre_quadrature_output_1

std::numeric_limits<local::float_type>::digits10: 101
laguerre_order: 2291

Finding the approximate roots...
root_estimates.size(): 1, 0.0%
root_estimates.size(): 8, 0.3%
root_estimates.size(): 16, 0.7%
...
root_estimates.size(): 2288, 99.9%
root_estimates.size(): 2291, 100.0%


Calculating abscissas and weights. Processed 1, 0.0%
Calculating abscissas and weights. Processed 9, 0.4%
...
Calculating abscissas and weights. Processed 2289, 99.9%
Calculating abscissas and weights. Processed 2291, 100.0%
//] [/gauss_laguerre_quadrature_output_1]

//[gauss_laguerre_quadrature_output_2

airy_ai_value  : 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874867334
airy_ai_control: 0.13529241631288141552414742351546630617494414298833070600910205475763353480226572366348710990874868323
airy_ai_value  : 0.11392302126009621102904231059693500086750049240884734708541630001378825889924647699516200868335286103
airy_ai_control: 0.1139230212600962110290423105969350008675004924088473470854163000137882588992464769951620086833528582
...
airy_ai_value
: 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331427451e-22
airy_ai_control: 3.8990420982303275013276114626640705170145070824317976771461533035231088620152288641360519429331426481e-22

Total... result_is_ok: true

//] [/gauss_laguerre_quadrature_output_2]


*/
